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SUMMARY 

An FFP algorithm has been developed based on the method of Lee et al* for the 
prediction of sound pressure level from low frequency high intensity sources. In order 
to permit accurate predictions at distances greater than 2km, new correction 
procedures have had to be included in the algorithm. Certain functions, whose 
Hankel transforms can be determined analytically, are subtracted from the depth 
dependent Green's function. The distance response is then obtained as the sum of 
these transforms and the FFT of the residual k dependent function. One procedure, 
which permits the elimination of most complex exponentials, has allowed significant 
changes in the structure of the FFP algorithm, which has resulted in a substantial 
reduction in computation time. 

1 . INTRODUCTION 

Sound pressure levels at large distances from a point source close to the ground have 
been predicted using ray based procedures 1 in enhanced zones and residue 
calculations 2 in strong shadow zones. In the published literature (see references 1 and 
2) these predictions have been shown to be approximately valid. The errors for the 
predictions in the enhanced zone increase when ground reflections become important 
and when landing ray densities become small. In the shadow zone errors increase 
when the sound speed gradient becomes small. Both the above procedures are 
inaccurate or inoperable in the transition regions between shadow and enhancement. 

Since the publication of the first paper on the FFP 3 for atmospheric sound 
propagation this method has been increasingly used for sound pressure level 
prediction. This has largely occurred because the FFP can operate irrespective of 
whether there are shadow or enhanced or even mixed conditions present. Moreover 
the FFP can take proper account of ground reflection. 

The most widely known FFP algorithm, the CERL-FFP, stems from the method of 
Lee et al 4 , which was a development of the algorithm described in Raspet et al's first 
paper 3 . 

Starting from reference 4 we have reworked the analysis to enable us to produce our 
own FFP algorithm, structured in such a way that we could incorporate a variety of 
corrections in k space and thereby extend the range of validity of the result in the 
transformed (real) space. 


* Lee et al. J. Acoust. Soc. Am., 79, 1986, pp 628-634. 
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In section 2 of this paper a brief description is given of the development of our first 
prototype algorithm. A survey of previously published k spectral corrections applied to 
the first prototype is given in section 3. Section 4 describes our second algorithm 
which is based on a '7 averaging* procedure and section 5 describes a technique 
adopted for speeding up the computation. In the last section, 7, a comparison 
between the FFP and our other model's predictions is given for a realistic case. 

2. DEVELOPMENT OF THE FIRST FFP ALGORITHM 

2.1 Sign Convention 

Lee et al 4 replaced the system of atmospheric strata by an analogue electrical 
network. In our reworking of the model we found that this was not necessary and 
that retention of acoustic equations for pressure and particle velocity ensured greater 
clarity. In addition our analysis showed that great care had to be exercised with the 
signs used in the ladder calculation. 

Raspet 3 correctly drew attention to the need for different signs for the characteristic 
admittance dependent on the direction of the particle velocity, which gives the correct 
sign for the characteristic admittance of top stratum, Y co = + iy 0 /o)p 0 where y g is 
the propagation constant and p Q is the density and for the characteristic admittance of 
the ground Y c m = + 1/Z where Z is the usual ground impedance. 

None of the published papers on the FFP, including the most recent ones (see Franke 
et al 5 ), make it clear that the signs used in the ladder admittance calculation must be 
chosen in accordance with the direction in which the calculation is performed. The 
equation for calculating an admittance (Y new ) at the stratum interface nearer to the 
source from the admittance at the stratum's other interface (Yqj^) is 

Yrm W « — — ? nh + Y Q l d / Y c m (1) 

'new 1 cm \ 1 / 

1 + € (^old/^cm) tan h Y m C m 


where the characteristic admittance for layer m, is Y cm = e iy m /o)p m , 7 m being the 
propagation constant, p m the density and £ m the stratum thickness; c is +1 when 
working upwards and -1 when working downwards. We note that the e's cancel in 
(1) so that the same equation can be used whether working above or below the 
source interface. However in view of e in the top semi-infinite layer being set to +1 
the admittance calculated just above the source interface, Yfz^ - , has opposite sign to 
that calculated just below that height, Y(z s ) + . 

The z dependent Green's function at the source, Pfz^, is obtained from the known 
discontinuity in particle velocity at that location. Using our sign convention and 
noting the opposite signs of Y^)" 1- and Y(z s )“ a negative sign must appear in the 
denominator of P(z s ). 


P(z s > 


- (2i/a>p g ) 

Y(z s )+ - Y(z s ) “ 


( 2 ) 


Likewise performing the calculation of the Green's function at the detector, P(zq)* 
and retaining the above convention for t, the equation for calculation of p at the 
stratum interface nearest the source (P new ) from that at that stratum's other interface 
(P Q id) becomes 


Pnew “ P 


old | cosh 7m C m + £ sinh 7m C m [ J J < 3 > 


irrespective of whether the receiver is above or below the source. 


2.2 Ground Impedance 

The early FFP algorithms used the Delany Bazley 6 model for the ground impedance. 
Attenborough has suggested that his four parameter model 7 be used instead since this 
gives much smaller and more realistic impedance values at the low frequencies of 
interest in this study. 

2.3 Damping Coefficient 

The k spectrum calculated with the above ladder procedure has a large spike at the k 
value closest to co/c Q where c 0 is the sound speed for the top semi-infinite top layer 
and also a number of subsidary spikes at k values nearest to co/c m . These infinities 
produce errors when the Fourier transform is performed. A global damping is usually 
introduced, preferably by subtracting a small imaginary quantity ip from k, where p is 
typically 10” 4 . The transform is corrected for the effect of the damping by 
multiplying it by e/ 0 * 

This procedure effectively produces different damping effects on each spike dependent 
on their proximity to the nearest k sample point. The k sampling errors are 
therefore only partially removed. 

3. APPLICATION OF PUBLISHED k SPECTRAL CORRECTIONS 


3.1 Candel and Crance's k Shift Procedure 

Candel and Crance 8 proposed a method which ensured the k sample intervals were 
chosen so that the main peaks all occurred well away from the sample points. 
Thereby the transform errors arising from the presence of spikes in the k spectrum 
were substantially reduced. 

We set up an algorithm to readjust the k sampling interval, Ak, until the 5 main 
spikes were more than a limiting distance (Ak/10) from the nearest k sample point. 
The method requires a full ladder calculation for typically l/8th of the total k range 
centred on lqj = oj/c 0 for each setting of Ak. This can be time consuming when 
large numbers of strata are considered. 

The method worked well for large Ak values when the total number of k samples, N, 
was less than lk. For the larger N values required for ranges greater than 2km it 
became increasingly difficult to ensure the main spikes were distanced more than the 
above limiting value from the sample points. 

* The term k spectrum in this paper refers to F = k P(zjy) which is the function 
to be Fourier transformed . 



3.2 The Richards and Attenborough 9 Tail Remover 

The k spectrum very often has a non zero asymptotic value or tail as k goes to 
infinity. This tail invariably occurs when source and receiver are approximately the 
same height above ground. 

Introduction of a cut-off for the spectrum at a value k max produces large ripples in 
the transform which increase with range and with the proximity of k max to gj/c q . 3 
In this study we set k max at a value where the change in F was less than a 
prescribed limit, but this remained unsatisfactory. 

The k spectral tail can be removed by subtracting the function g(k) given by Richards 
and Attenborough 9 from F(k). 


g(k) - (A + e kz s) (1 - e ~ 6k ) 


(4) 


where Zg is the source height and A = \F (k max ) | and 5 is the derivative of 
\F (k)i at small k. g(k) has an exact Hankel transform which is added to the 
Fourier transform of F (k) - g (k). The above tail remover does not have the 
correct phase at small k and this produces a small error in the transform. 

Experiments with functions, with exact Hankel transforms which mimic the F (k) 
behaviour at small k and near the main spike are in progress. 

In Figure 1 the attenuation for a source and receiver 2m above an impedance ground 
at 50Hz is shown with and without the tail remover. 

4- MODIFIED FFP ALGORITHM USING 'y' AVERAGING 1 

For the large ranges of interest in this study very small Ak values must be used 
therefore the Candel and Crance method does not work. A novel technique for 
obtaining a k spectrum which can be more reliably transformed has been developed. 

The y averaging procedure is most easily understood by considering a single step in 
the admittance calculation for one stratum, index m, as described in section 2. 

Writing equation (1) in matrix form 


P 
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new 


P 


U 
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(5) 


where the matrix M m has elements which are functions of k and are usually evaluated 
at a single k sample value, say k r . 


M m- 


esinh y m Q m 
^cm 


cosh 7m C m 


( 6 ) 


where 7m -f k r 9 - u>Vc m 2 
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For Fourier integrals of smoothly varying functions F(k) (which have to be 
approximated by a discrete Fourier transform), the optimal sampling is at equidistant 
k r values. These F(k) samples can be assumed to approximate the contribution to the 
integral over the interval k r - Ak/2 to k r + /ik/2. In the presence of integrable 
infinities a different averaging process must be used. 

The simplest method follows from changing the integration variable in the above 
interval from k to 7; then elementary averaging leads to (6) with 7 m (k r ) replaced 

by 

Tm ” [ 7m O^r ~ + Tm Ohr + Ar/2) J/2 

We employ this averaging method for those intervals where 7 2 changes sign. 

Provided c m is real (no damping) the values of y m in the integrand are either pure 
imaginary or pure real and so are the 7 m values at all sample points except those 
where the above averaging is employed. Hence the hyperbolic functions in (6) can be 
replaced by real trigonometric or hyperbolic functions, the full complex functions being 
required at the above critical points only. A speed up by a factor of 3 was obtained 
by this procedure over that using (6) in its general complex form. This high speed 
algorithm cannot be used if artificial damping is present. 

5. k SPACE INTERPOLATION 

In order to obtain predictions out to large distances the number of k samples, N, for 
a given k,jj ax must be increased. As the major part of the computation lies in the 
determination of the z dependent kernel at each k value, we can achieve considerable 
improvements in calculation speed by interpolating F(k) where it varies smoothly. 

This applies to all regions of the spectrum lying outside the range spanned by oi/c m 
for all layers (typically k max /8). It was found that F(k) only needed to be evaluated 
every 8th point in the smooth region. 

For most cases linear interpolation of F(k) is adequate. However at very large ranges 
fluctuations occur in the transform due to the small discontinuities in the slope of the 
interpolated F(k). These errors can be reduced by using cubic interpolation. 

6. STRATUM QUANTISATION ERROR 

In Figure 2 the predicted attenuation above an impedance ground is shown for two 
different stratum configurations with the same small linear sound speed gradient (0.01 
s -1 ). For the solid curve all strata above 30m are taken as 27m thick and for the 
dotted curve as 54m thick. It is clear that the undulations in the dotted curve are 

much bigger than for the solid curve. This occurs because the FFP uses mid stratum 

sound speed averages for the whole of each stratum. This produces a stratum 

sampling error which gets worse the thicker the strata and the larger the sound speed 

gradient. There are two remedies: one is to use many very fine strata that increases 

the amount of computation, the other is to use linear sound speed variations in each 
stratum and use a modified procedure employing Airy functions 1 °- 

7. COMPARISON OF THE FFP PREDICTIONS WITH 

RAY/RESIDUE MODEL PREDICTIONS 

In Figure 3 the attenuation predicted by the FFP for a lapse-inversion-lapse 
meteorology is compared with that obtained by a hybrid method. This latter method 
combines the predictions from a ray based model for the enhanced regions with 
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those from a residue model for the shadow regions. The advantages of using the 
FFP are apparent in this figure. The residue model overpredicts the shadow 
attenuation and the ray model is restricted to giving predictions only in the region 
where there are ray landing points. Neither of these models can properly deal with 
the transitional region between shadow and enhancement. 
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Comparison with Hybrid Model Predictions 
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Figure 3 


Comparison of FFP Predictions with those from a Hybrid Model. 
Sound speed gradient up to 100m -0.1, from 100m to 200m +0.1 
and above 200m -0.05 and an impedance ground. 
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